Introduction

In natural environments or agricultural practices, growth and production of plants were usually suppressed by environmental stresses. Freeze injury caused by subzero temperature is one of the most common stresses, especially for plants distributed in sub-tropical regions (Rooy et al. 2017). When plants suffered from low temperature, various physiological processes can be activated, including molecular and biochemical processes, and by altering gene expression level and phenotype to adapt to low temperature. In previous studies, low temperature down-regulated photosynthesis-related protein expression, sugar transport, respiratory rate, while up-regulated carbohydrate metabolism, cell wall remodeling, redox adjustment, and defence/detoxification (Janmohammadi et al. 2015; Ding et al. 2019; Leites et al. 2019). Abscisic acid (ABA) plays a very crucial role for improving plant cold resistance (Fu et al. 2017). Protein biosynthesis induced by low temperature was similar with ABA or light water shortage, and these proteins possessed the functions of anti-freezing and osmotic adjustment (Mantyla et al. 1995). It is now possible to investigate the response of plant to environmental stresses on global transcriptional level with the development of transcriptomics.

Phoebe zhennan, commonly called Zhennan, is an evergreen plant that belongs to the Lauraceae family and is mainly distributed in southwest of China (Chen et al. 2018). Zhennan is a high-quality timber tree species and is well known for its remarkable surface golden tint. Since Ming Dynasty, the wood was widely used as the materials of palace buildings, such as the Forbidden City (Xie et al. 2015). This tree species is also characterized as an ornamental plant, usually for landscaping (Hu et al. 2015). However, the extensive deforestation in recent centuries, the large diameter Zhennan timbers are increasingly scarce. In the past few years, the natural resources of Zhennan have been continuously declined (Cao et al. 2016; Gao et al. 2016). In the recent years, the conservation and recovering of rare and precious tree species resources, including Zhennan, received more attention (Hu et al. 2015). At present, cultivation of Zhannan has gained attention for ecological restoration and use in timber resource. To our knowledge, introduction and cultivation of Zhennan has been carried out in southeast regions of China in the past few decades. As a sub-tropical tree species, however, Zhennan usually meet with freeze injury in winter when temperature below freezing.

Calcium ions (Ca²⁺), the crucial second messengers, involve in mediating responses to many abiotic stimuli, including low temperature, salinity and drought (Perochon et al. 2011; Wang et al. 2016). Previous studies found that the stress tolerance of plants can be improved through the Ca²⁺ application. Huang et al. (2015) found that exogenous Ca²⁺application enhanced the cold tolerance of tea plant. Trapnell et al. (2011) also found a post-drought recovery phenomenon of tea plant by foliar-sprayed with Ca²⁺. However, little information is available about the effects of exogenous Ca²⁺application on cold tolerance of Zhennan. The present study aims at investigating the physiological mechanisms of cold tolerance in Zhennan on transcriptional and metabolic levels. The information would be of great value for breeding and cultivation.

Materials and Methods

Plant materials and experimental design

The 1-year-old seedlings of Zhennan were cultivated in a controlled phytotron with a 50% relative humidity, 12 h photoperiod, and a 17°C/12°C diurnal/night temperature for 25 d. The leaves of Zhennan seedlings were then continuously sprayed daily with 10 mmol/L CaCl₂solution or distilled water for 7 d. On the eighth day, all seedlings were moved to open field and exposed to natural low temperature condition. The first leaves below the top buds were collected, and immediately frozen in liquid nitrogen, and stored at -70℃ before and after the cold treatment (i.e., 3 and 15 d). During the experimental period, the daily air temperature was shown in Fig. 1. There were 6 seedlings in treatments of CaCl₂-pretreated and control and three biological replicates in each treatment.

Measurements of physiological and biochemical parameters

The leaves of Zhannan seedlings were sampled according to each replication and the fresh leaves were used to measure physiological and biochemical parameters. The levels of ABA, peroxidase (POD), catalase (CAT) were measured based on the protocol of Plant ABA ELISA Kit, Plant POD ELISA Kit and Plant CAT ELISA Kit, respectively. The contents of chlorophyll were detected by using spectrophotometry at 663 and 645 nm, respectively, after extracting with acetone-ethyl alcohol (v/v = 17/3) solution (Paul and Driscoll 1997). Relative electrical conductivity was analyzed as a method described by Cui et al. (2018). The contents of proline, malondialdehyde (MDA), and soluble sugar were measured according to a method described by Bates et al. (1973), Draper and Hadley (1990), Buysse and Merckx (1993), respectively.

RNA extraction, cDNA library construction, and RNA-sequencing

In order to perform the transcriptome sequencing analysis, the leaves of Zhennan from CaCl₂-pretreated and control were sampled after 3 d-exposing to low temperature. A Picopure^TM RNA isolation Kit (Thermo Fisher Scientific, Waltham, U.S.) was used to extract the total RNA from the leaves of CaCl₂-pretreated and control of Zhennan, and then followed by RNA purification with RNase-free DNase I (TaKaRa, Dalian, China) and the tests of RNA degradation and contamination by using a 2% agarose gel electrophoresis. In addition, RNA concentration, purity and integrity were measured by using a Qubit 2.0 Fluorometer (Life Technologies, C.A., U.S.A.), a NanoPhotometer spectrophotometer (Implen, C.A., U.S.A.) and an Agilent 2100 analyzer (Agilent Technologies, C.A., U.S.A.), respectively. A modified procedure described by Yang and Huang (2018) was used to construct the cDNA libraries. Transcriptome sequencing of cDNA libraries was performed on an Illumina HiSeq 2500 platform, and 125 bp paired-end reads were generated.

Transcriptome assembly and functional unigene annotation

The Trinity de novo assembly program was used to assemble the remaining high-quality sequencing clean data into contigs after filter the low-quality reads (Feng et al. 2016). The BLAST software was then used to functionally annotate all assembled unigenes to the available sequences in public databases, including eggNOG4.5 (Huerta-Cepas et al. 2015), RefSeq non-redundant proteins (NR) (Deng et al. 2006), euKaryotic Orthologous Groups (KOG) (Koonin et al. 2004), Swiss-Prot (Apweiler et al. 2004), Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa et al. 2004), Gene Ontology (GO) (Ashburner et al. 2000) and Clusters of Orthologous Groups (COG) (Tatusov et al. 2000). The KEGG orthology was obtained by using the KOBAS 2.0 program, and then then the amino acid sequence of all unigenes were predicted. In addition, HMMER software was used to align to Protein family (Pfam) database for obtaining the annotation information of all unigenes (Finn et al. 2004).

Differential expression and enrichment analysis

After mapping the clean data back onto the assembled transcriptome with TopHat software, the expression levels of unigenes were calculated based on the fragments per kilobase of transcript per million mapped reads (FPKM) by using Expectation maximization (RSEM) (Li and Dewey 2011; Langmead and Salzberg 2012). And then, the differentially expressed genes (DEGs) between the control and the CaCl₂-pretreated seedlings of Zhennan were selected by using the DESeq R package (Anders and Huber 2010). In this study, the DEGs were characterized as the genes with an absolute log₂(fold change) value ≥ 2 and a false discovery rate (FDR) ≦ 0.001. Based on the Wallenius’ noncentral hypergeometric distribution, the GO enrichment analysis of DEGs was conducted with GOseq R package (Young et al. 2010). Furthermore, the KEGG pathway enrichment analysis of the DEGs was performed by using KOBAS software for discerning the related biochemical and signal transduction pathways (Mao et al. 2005).

qRT-PCR validation

The sequencing data was validated with a quantitative real-time polymerase chain reaction (qRT-PCR) analysis. In this study, 6 candidate genes were selected, and the gene-specific primers of these genes were designed by using the Primers Premier 5.0 program (Suppl. Table S1). After synthesizing cDNA by using the above isolated total RNA, the qRT-PCR was performed on a StepOnePlus Real-Time PCR Systems (Applied Biosystems, CA, USA) with the following PCR procedure conditions: 95℃ for 120 s, 40 cycles of 95℃ for 10 s, 60℃ for 30 s, 72℃ for 30 s. The relative expression levels of the selected genes normalized against the two Zhennan gene (pentatricopeptide repeat-containing protein At1g43980 and DNA polymerase epsilon catalytic subunit A-like) expression leaves and were calculated by using the 2^-^ΔΔ^C_T method (Livak and Schmittgen 2001). There were 3 biological replications for each selected gene and 3 technical replications for each biological replicate.

Statistical analysis

S.P.S.S. version 16.0 (S.P.S.S. Inc., Chicago, I.L., U.S.A.) was used to perform the analysis of variance (ANOVA) and the significant differences among the treatments were calculated with Duncan's multiple-range test. All statistical analyses were performed at a 95% confidence level.

Results

Effects of exogenous CaCl₂ on physiological and biochemical characteristics in cold-stressed Zhennan leaves

Eight physiological parameters were analyzed to clarify the effects of exogenous CaCl₂on Zhennan cold resistance. CaCl₂-pretreatment significantly increased the ABA content in Zhennan leaves (P < 0.05), while the control group possessed a significant higher ABA content after exposed to natural low temperature (Fig. 2a). The chlorophyll content, however, had an adverse variation pattern during the experimental period compared with ABA in leaves (Fig. 2b). Activities of the two peroxidases (POD and CAT) were significantly enhanced by exogenous CaCl₂ before Zhennan seedlings were exposed to low temperature (0 d, P < 0.05, Fig. 2c–d). After Seedlings were exposed to low temperature environment (3 and 15 d), the CaCl₂ treatment still had a relative higher peroxidase activity, though no statistically significant was detected. Throughout the experiment, relative electrical and the contents of MDA, soluble sugar, and proline were significantly increased (P < 0.05, Fig. 2e–h). Seedlings in control group had a significant higher MDA content and relative electrical, while CaCl₂-treated seedlings had a significant higher proline and soluble sugar contents. Furthermore, the shoot tips of control seedlings showed mild

Fig. 1: The daily mean temperature and daily minimum temperature during the experimental periods (from December 6 to 21, 2018)

Fig. 2: Effects of exogenous CaCl₂ pretreatment on abscisic acid (ABA) (a), chlorophyll (b), peroxidase (POD) (c), catalase (CAT) (d), malondialdehyde (MDA) (e), relative electrical conductivity (f), proline (g), and soluble sugar (h) in cold-stressed Zhennan leaves. Seedlings of Zhennan were sprayed with distilled water (control) and 10 mmol/L CaCl₂solution for 7 d under normal growth condition. Values within each graph followed by the different letters indicate significant differences between treatments (lower case), and among different sampling time (upper case) (n = 3) according to Duncan’s test (P < 0.05)

wilt, while CaCl₂-treated seedling displayed normal growth (Fig. 3). These data indicated that CaCl₂ pretreatment increased the cold tolerance of Zhennan seedlings.

Table 1: Summary of transcriptome sequencing data

Length range	Transcript	Unigene
200-300	36,992 (18.32%)	30,857 (30.79%)
300-500	27,238 (13.49%)	15,350 (15.32%)
500-1000	45,082 (22.32%)	19,673 (19.63%)
1000-2000	51,015 (25.26%)	18,768 (18.73%)
≥ 2000	41,643 (20.62%)	15,555 (15.52%)
Total number	201,970	100,203
Total length	257,180,926	104,144,175
N50 length	2,035	1,910
Mean length	1273.4	1039.3

Fig. 3: Phenotypes of Zhennan seedlings under different growth conditions throughout the 15 d low temperature treatment

Transcriptome sequencing, assembly and functional annotation

To clarify the molecular regulation mechanism of CaCl₂ on cold resistance of Zhennan, transcriptome sequencing was performed on a Illumina Hiseq instrument system. After filtering out the low-quality data, approximately 40 726 813 and 54 928 919 clean reads were obtained for the 3 d-CaCl₂ and 3 d-control, respectively (Suppl. Table S1). Data quality assessment showed that > 94.33% of these data possessed a quality score of Q30 (i.e., the sequencing error rates < 0.1%), and the GC contents for CaCl₂-pretreated and control were 46.18%, and 47.48%, respectively (Suppl. Table S1). Furthermore, 29 166 060 (71.61%) and 39 352 390 (71.64%) clean reads in two cDNA libraries that could be successfully matched to the reference genome. The Zhennan transcriptome was assembled by using Trinity program, the final assembled 100 203 unigenes had a mean length of 1039.3 bp and an N50 length of 1910 bp (Table 1). These results indicated that the unigenes were qualified for further functional annotation.

The unigenes were functionally annotated with eight public databases (GO, Swiss-prot, NR, KOG, COG, KEGG, eggNOG and Pfam) for analyzing the gene function information. The unigenes of Zhennan were annotated as a transcriptome and a total of 47 781 unigenes were matched in NR database, in which 47 716 unigenes were homologous to multiple species genomes, including Nelumbo nucifera (11.99%), Macleaya cordata (7.01%), Vitis vinifera (3.76%), Quercus suber (3.43%), Hortaea werneckii (2.70%), Elaeis guineensis (2.66%), Alternaria alternate (2.36%), Phoenix dactylifera (2.22%), Baudoinia panamericana (2.04%), Aquilegia coerulea (1.94%), and other (59.90%) (Suppl. Fig. S1).

For classifying genes and gene products, 20 712 unigenes were divided into 3 GO categories and 50 functional subcategories (Fig. 4 and Suppl. Table S2). For the category of cellular component, the top 6 subcategories were cell (18039 unigenes, 87.1%), membrane (17149 unigenes, 82.8%), organelle (16703 unigenes, 80.6%), membrane part (15991 unigenes, 77.2%), organelle part (14476 unigenes, 69.9%), and macromolecular complex (14031 unigenes, 67.7%). For the molecular function category, the most highly represented 4 GO terms were catalytic activity (18436 unigenes, 89.0%), binding (18208 unigenes, 87.9%), transporter activity (12518 unigenes, 60.4%) and structural molecule activity (10925 unigenes, 52.7%). For biological process, the dominant groups were metabolic process (18459 unigenes, 89.1%), cellular process (18231 unigenes, 88.0%), single-organism process (17116 unigenes, 82.6%), response to stimulus (13656 unigenes, 65.9%), localization (14544 unigenes, 70.2%), biological regulation (14680 unigenes, 70.9%) and cellular component organization or biogenesis (13429 unigenes, 64.8%). Our results indicated that membrane component-related genes were assigned to the cellular component category (GO: 0016021), while the genes related to oxidoreductase activity (GO: 0016491) and positive regulation of abscisic acid-activated signaling pathway (GO: 0009789) were enriched to biological process.

The COG database was also used to evaluate the homology of gene products. A total of 14 808 unigenes were classified into 26 COG categories (Fig. 5 and Suppl. Table S3). Carbohydrate transport and metabolism (G, 10.01%), chaperones, posttranslational modification, protein turnover (O, 9.73%), metabolism and lipid transport (I, 7.81%), metabolism and amino acid transport (E, 7.64%), signal transduction mechanism (T, 7%), and cell wall/membrane/envelope biogenesis (M, 4.03%) contained 1469, 1428, 1146, 1122, 1027, and 592 unigenes, respectively. This result indicated that these unigenes might be involved in cold tolerance of Zhennan seedlings.

Identification and analysis of DEGs

According to the DESeq program, 4245 unigenes were characterized as the DEGs between the paired samples. Within these DEGs, 477 unigenes were up-regulated, while 3768 unigenes were down-regulated for the samples exposed to low temperature compared with control group (Suppl. Fig. S2). These DEGs were functionally annotated into the 8 public databases to further analyse the function of the DEGs: 3776 in NR, 3600 in eggNOG, 3101 in Pfam, 2548 in GO, 2342 in KOG, 2105 in Swiss-Prot, 1680 in COG, and 1628 in KEGG, respectively. Based on GO functional annotation, 2548 DEGs were assigned to three categories: molecular function, cellular component and biological process (Fig. 4). Among the category of biological process, the top five subcategories were metabolic process (2308 DEGs), cellular process (2285 DEGs), single-organism process (2157 DEGs), localization (1946 DEGs), and biological regulation (1894 DEGs). Within cellular component category, cell (2305 DEGs), cell part (2304 DEGs), membrane (2301 DEGs), organelle (2199 DEGs), and membrane part (2091 DEGs) were the most highly represented. Catalytic activity, binding, and transporter activity were the dominant subcategories in molecular function group.

Fig. 4: Classification of all assembled unigenes and DEGs annotated in GO database. This figure shows the ontology enrichment of DEGs in various GO functional subcategories. All DEGs fell into the 3 categories of biological process, celluar component and molecular function.

Fig. 5: Classification of all assembled unigenes in COG database. The horizontal axis refers to the 26 categories, and the vertical axis indicates the number of DEGs in each category

The KEGG database was used for the systematic analysis of gene function. In the present study, a total of 17 376 unigenes and 1766 DEGs were classified into 121 known pathways. The top 50 KEGG pathways were divided into four categories: 3 pathways related to cellular processes, 1 pathway connected with environmental information processing, 13 relevant to genetic information processing, and 33 pathways associated to metabolism (Fig. 6).

Validation of DEGs by qRT-PCR

To further verify the reliability of the RNA-Seq data, six cold-responsive DEGs were randomly selected and analyzed with qRT-PCR (Suppl. Table S3). The candidate genes Text Box:

Fig. 6: The KEGG pathway assignment of the DEGs. The vertical axis refers to the KEGG metabolic pathways, and the horizontal axis indicates the number of DEGs in each pathway encode RPR (c183911), Rboh (c167842), RPM1 (190929), RPS2 (173175), HSP90 (180553) and PR1 (c192795), respectively. As shown in Fig. 7, the log₂(fold change) values obtained by qRT-PCR were closely consistent with the transcript abundance detected by RNA-Seq for all of the selected unigenes, despite the difference in the absolute fold change between the two methods. This result confirmed that the reliability of the RNA-Seq analysis.

Discussion

Plants have evolved elaborate mechanisms of resistance to cold stress (Ding et al. 2019). At the biochemical and physiological levels, numerous protective substances are biosynthesized in plants, such as cold-resistance proteins, proline, and soluble sugar (Kaplan et al. 2007). These substances are involved in ice crystal formation, reactive oxygen species (ROS) scavenging, osmotic potential, and regulating the stability of cell membranes (Dong et al. 2009). Of these substances, ABA plays a very important role in frost resistance induction of plants. In the previous studies, plants acquire increased cold tolerance upon prior drought acclimation, a process by which the endogenous ABA in leaves increased (Mantyla et al. 1995; Gusta et al. 1996). In addition, exogenous ABA treatment can increase cold tolerance of plants under a non-cold acclimation temperature condition (Gusta et al. 1996). In this study, we found a significant higher ABA content in leaves of Zhennan before seedlings were exposed to low temperature, which indicated that CaCl₂ pretreatment has a similar effect to cold acclimation. The two low-molecular-weight solutes, proline and soluble sugars, function as osmotic adjustment to protect plants against the damage caused by cold stress (Ruelland et al. 2009). While MDA and relative electrical conductivity are usually used to reflect lipid peroxidation, and the elevated value suggests the membrane systems may be severely damaged (Catola et al. 2016). In the present study, the increased proline and soluble sugars, while decreased MDA and relative electrical conductivity indicated that CaCl₂-treated seedlings of Zhennan had a higher cold resistance, which was consistent with the previous report (Tan et al. 2011). The CaCl₂ pretreatment also enhanced the activities of peroxidases (POD and CAT) before seedlings were exposed to low temperature (Fig. 2c–d). This result agreed with the lower MDA and relative electrical conductivity for CaCl₂-treated Zhennan.

The cold damaged plants usually showed photosynthetic inhibition, blocked protein synthesis, and desaturation of membrane lipid (Williams et al. 1988; Shinozaki and Yamaguchi-Shinozaki 2000; Rajkowitsch et al. 2007). Transcriptome sequencing was performed to clarify the molecular regulation mechanism of CaCl₂ on Zhennan cold resistance. To our knowledge, none of the literatures have dealt with the transcriptome sequencing analysis on Zhennan. In this study, DEGs analysis revealed that abundant biological processes associated to cold resistance have been changed by low temperature treatment. For example, calcium-transporting ATPase activity (GO: 0005388), Ca²⁺ transmembrane transport (GO: 0070588), peroxidase activity (GO: 0004601), response to oxidative stress (GO: 0006979), negative regulation of ABA-activated signaling pathway (GO: 0009788) and cell surface receptor signaling pathway (GO: 0007166) were markedly down-regulated. While, calmodulin binding (GO: 0005516), protein serine/threonine kinase activity (GO: 0004674) and protein phosphorylation (GO: 0006468) were significantly upregulated.

For organisms, biological function links closely with the coordination of different genes. It is very important to further reveal the gene function through analyzing the pathway of DEGs. The pathways of phagosome (18 DEGs), peroxisome (25 DEGs), and endocytosis (43 DEGs) are associated to antioxidant defense (Liu et al. 2018). The phosphatidylinositol signaling system is related to signal transduction. About 78.3% DEGs were assigned to metabolism, including fatty acid metabolism, amino acid metabolism, carbon and carbohydrate metabolism. Furthermore, several pathways are involved in secondary metabolism, such as biosynthesis of flavonoid, terpenoid, and alkaloid. Overall, DEGs annotated into above pathways play important role for Zhennan seedlings resisting cold stress.

The systematic defense responses to biotic and abiotic stresses are the important resistance mechanism for plants (Bonello et al. 2006). In the present study, the low temperature activated the plant-pathogen interaction pathway (ko04626) in Zhennan seedlings (Fig. 8). The previous study has found several cold-induced antifreeze proteins in leaves of cold-acclimation plant (Thomashow 2001). The antifreeze proteins were also found can be induced by different pathogen infection, and they are known as pathogenesis-related protein (PR). In plants, these antifreeze proteins possess the dual function of cold and disease resistances, so as to protect plants from cold stress and pathogen infection. In the present study, 11 DEGs were found in the plant-pathogen interaction pathway based on the KEGG enrichment analysis. Among them, respiratory burst oxidase homolog protein (Rboh, c167842), calcium-dependent protein kinase (CDPK, c173040) and disease resistance protein RPM1 (c190929) were up-regulated, while pathogenesis-related protein 1 (PR1, c173293, c192795), disease resistance protein RPS2 (c173175), and molecular chaperone HtpG (HSP90, c180553) were down-regulated. In addition, genes of calmodulin/calcium binding protein CML (CaM/CML, c164829, c191287/c162995, c184878) showed mixed regulation.

As a second messenger, calcium ion is very important for plant to response to environmental stimuli. The rapidly enhanced cytosolic Ca²⁺ concentration via Ca²⁺ channels after cold stress is considered as one of the earliest cold signaling events (Anders and Huber 2010).Interestingly, Knight et al. (1996) found that cold-induced COR gene was Ca²⁺-dependent. The more recently study indicated that G-protein regulator involved in cold sensing through modulating calcium signals in rice (Oryza sativa) (Ma et al. 2015). These studies reveal that Ca²⁺plays a key role in various biological processes, including cold stress sensing, signal transduction, and regulation of gene expression. In addition, numerous evidences indicate that there is a strong linkage between Ca²⁺ and other messenger molecules such as nitric oxide (NO) and ROS. For example, low ROS concentration can promote Ca²⁺ influx into the cytoplasm. On the other hand, Ca²⁺ can regulate ROS production in plant under various stresses, such as low temperature, heat, and drought stress (Wang et al. 2016). As such, exogenous Ca²⁺ application may enhance the cold resistance of plants through activating the ion channels in the cell membrane and further inducing

Fig. 8: The KEGG pathway of the plant-pathogen interaction. The red boxes represent these enzymes are related to up-regulated DEGs, the green boxes refer to the down-regulated DEGs, and the blue box indicates the enzymes were both associated to up- and down-regulated DEGs

Fig. 7: Validation of the relative expression levels of the 6 randomly selected DEGs by qRT-PCR. Data refer to the log₂(fold change) of expression for each selected DEG in the cold treatment compared to control group

expression of resistance genes. Indeed, in the present study, the cold tolerance of Zhennan seedlings were improved by exogenous Ca²⁺ pretreatment, and this phenomenon have been proved on the levels of phenotype and physiological characteristics (Fig. 2 and 3). Also, the gene expression in many metabolic pathways, including plant-pathogen interaction pathway, have been changed by exogenous Ca²⁺ pretreatment.

Compared with control seedling, exogenous Ca²⁺ pretreatment strongly activated the genes of PzCDPK and PzRboh, resulting in the accumulation of ROS by PAMP-triggered immunity and finally the defense responses were activated by hypersensitive response (HR) and cell wall reinforcement. The mixed regulated genes of PzCaM/CML and disease resistance protein (PzRPM1, PzRPS2) were not resulting in the change of downstream genes. Interestingly, the defense-related genes induced by WRKY transcription factor families (WRKY TFs) were not change or down-regulated. In other plant species, WRKY TFs have been widely reported that they can enhance the resistance of plant to various stresses and pathogen infection (Haider et al. 2017). These results indicated that CDPK might play key role for improving the cold resistance of Zhennan seedlings.

Conclusion

In conclusion, our results indicated that the exogenous CaCl₂-pretreatment increased the cold resistance of P. Zhennan seedlings, which was verified by phenotype and physiological characteristics. The following comparative transcriptome analysis revealed that numerous genes in cold resistance-related biological processes were activated, such as calcium-transporting ATPase activity, peroxidase activity, negative regulation of ABA-activated signaling pathway, etc. In addition, the KEGG enrichment analysis suggested that the CDPK within plant-pathogen interaction pathway might play key role for improving the cold resistance of P. zhennan. The findings provide a novel insight into the complexity of the Ca²⁺-mediated cold resistance of P. Zhennan seedlings and can improve cold stress tolerance of this plant.

Acknowledgements

We acknowledge financial support from the National Key Research and Development Program of China (2016YFD0600603).

Author Contributions

Bo Deng performed measurements, analyzed data and wrote the manuscript. Botao Jia conducted the experiment, performed measurements and analyzed data. Guihua Liu conceived the experimental idea and provided field support in establishing the experiment. Xiaoyan Zhang established the experiment and performed measurements.

References

Anders S, W Huber (2010). Differential expression analysis for sequence count data. Genome Biol 11; Article R106

Apweiler R, A Bairoch, CH Wu, WC Barker, B Boeckmann, S Ferro, E Gasteiger, HZ Huang, R Lopez, M Magrane, MJ Martin, DA Natale, CO Donovan, N Redaschi, LL Yeh (2004). UniProt: The universal protein knowledge base. Nucl Acids Res 32:115‒119

Ashburner M, CA Ball, JA Blake, D Botstein, H Butler, JM Cherry, AP Davis, K Dolinski, SS Dwight, JT Eppig, MA Harris, DP Hill, L Issel-Tarver, A Kasarakis, S Lewis, JC Matese, JE Richardson, M Ringwald, GM Rubin, G Sherlock (2000). Gene ontology: Tool for the unification of biology. Nat Genet 25:25‒29

Bates LS, RP Waldren, ID Teare (1973). Rapid determination of free proline for water-stress studies. Plant Soil 39:205‒207

Bonello P, TR Gordon, DA Herms, DL Wood, N Erbilgin (2006). Nature and ecological implications of pathogen-induced systemic resistance in conifers: A novel hypothesis. Physiol Mol Plant Pathol 68:95‒104

Buysse J, R Merckx (1993). An improved colorimetric method to quantify sugar content of plant tissue. J Exp Bot 44:1627‒1629

Cao J, H Shang, Z Chen, Y Tian, H Yu (2016). Effects of elevated ozone on stoichiometry and nutrient pools of Phoebe bournei (Hemsl.) Yang and Phoebe zhennan SLee et FNWei seedlings in subtropical China. Forests 7; Article 78

Catola S, G Marino, G Emiliani, T Huseynova, M Musayev, Z Akparov, BE Maserti (2016). Physiological and metabolomic analysis of Punica granatum (L.) under drought stress. Planta 243:441‒449

Chen Z, J Cao, H Yu, H Shang (2018). Effects of elevated ozone levels on photosynthesis biomass and non-structural carbohydrates of Phoebe bournei and Phoebe zhennan in subtropical China. Front Plant Sci 9; Article 1764

Cui J, N Jiang, XX Zhou, XX Hou, GL Yang, J Meng, YS Luan (2018). Tomato MYB49 enhances resistance to Phytophthora infestans and tolerance to water deficit and salt stress. Planta 248:1487‒1503

Deng YY, JQ Li, SF Wu, YP Zhu, YW Chen, FC He (2006). Integrated nr database in protein annotation system and its localization. Comput Eng 32:71‒74

Ding YL, YT Shi, SH Yang (2019). Advance and challenges in uncovering cold tolerance regulatory mechanisms in plants. New Phytol 222:1690-1704

Dong CH, BK Zolman, B Bartel, BH Lee, B Stevenson, M Agarwal, JK Zhu (2009). Disruption of Arabidopsis CHY1 reveals an important role of metabolic status in plant cold stress signaling. Mol Plant 2:59‒72

Draper HH, M Hadley (1990). Malondialdehyde determination as index of lipid-peroxidation. Meth Enzymol 186:421‒431

Feng YZ, L Zhang, JM Fu, FD Li, L Wang, XF Tan, WJ Mo, HP Cao (2016). Characterization of glycolytic pathway genes using RNA-Seq in developing kernels of Eucommia ulmoides. J Agric Food Chem 64:3712‒3731

Finn RD, A Bateman, J Clements, P Coqqill, RY Eberhardt, SR Eddy, A Heger, K Hetherington, L Holm, J Mistry, EL Sonnhammer, J Tate, M Punta (2004). Pfam: The protein families database. Nucl Acids Res 42:222‒230

Fu JJ, Y Wu, YJ Miao, YM Xu, EH Zhao, J Wang, HE Sun, Q Liu, YW Xue, YF Xu, TM Hu (2017). Improved cold tolerance in Elymus nutans by exogenous application of melatonin may involve ABA-dependent and ABA-independent pathways. Sci Rep 7; Article 39865

Gao J, W Zhang, J Li, H Long, W He, X Li (2016). Amplified fragment length polymorphism analysis of the population structure and genetic diversity of Phoebe zhennan (Lauraceae) a native species to China. Biochem Syst Ecol 64:149‒155

Gusta LV, RW Wilen, P Fu (1996). Low-temperature stress tolerance: The role of abscisic acid sugars and heat-stable proteins. Hortic Sci 31:39‒46

Haider MS, MM Kurjogi, M Khalil-Ur-Rehman, M Fiaz, T Pervaiz, ST Jiu, HF Jia, W Chen, JG Fang (2017). Grapevine immune signaling network in response to drought stress as revealed by transcriptomic analysis. Plant Physiol Biochem 121:187‒195

Hu Y, B Wang, TX Hu, H Chen, H Li, W Zhang, Y Zhong, HL Hu (2015). Combined action of an antioxidant defence system and osmolytes on drought tolerance and post-drought recovery of Phoebe zhennan SLee saplings. Acta Physiol Plantarum 37; Article 84

Huang YT, WJ Qian, B Wang, HL Cao, L Wang, XY Hao, XC Wang, YJ Yang (2015). Effects of exogenous calcium and inhibitors of calcium signaling transduction pathway on cold resistance of tea plant. J Tea Sci 35:520‒526

Huerta-Cepas J, D Szklarczyk, K Forslund, H Cook, D Heller, MC Walter, T Rattei, DR Mende, S Sunagawa, MKLJ Jensen, CV Mering, P Bork (2015). eggNOG 4.5: A hierarchical orthology framework with improved functional annotations for eukaryotic prokaryotic and viral sequences. Nucl Acids Res 44:286‒293

Janmohammadi M, L Zolla, S Rinalducci (2015). Low temperature tolerance in plants: Changes at the protein level. Phytochemistry 117:76‒89

Kanehisa M, S Goto, S Kawashima, Y Okuno, M Hattori (2004). The KEGG resource for deciphering the genome. Nucl Acids Res 32:277‒280

Kaplan F, J Kopka, DY Sung, W Zhao, M Popp, R Porat, CL Guy (2007). Transcript and metabolite profiling during cold acclimation of Arabidopsis reveals an intricate relationship of cold-regulated gene expression with modifications in metabolite content. Plant J 50:967‒981

Knight H, AJ Trewavas, MR Knight (1996). Cold calcium signaling in Arabidopsis involves two cellular pools and a change in calcium signature after acclimation. Plant Cell 8:489‒503

Koonin EV, ND Fedorova, JD Jackson, AR Jacobs, DM Krylov, KS Makarova, R Mazumder, SL Mekhedov, AN Nikolskaya, BS Rao, LB Rogozin, S Smirnov, AV Sorokin, AV Sverdlov, S Vasudevan, YI Wolf, JJ Yin, DA Natale (2004). A comprehensive evolutionary classification of proteins encoded in complete eukaryotic genomes. Genome Biol 5:1-28

Langmead B, SL Salzberg (2012). Fast gapped-read alignment with Bowtie 2. Nat Meth 9:357‒359

Leites LP, GE Rehfeldt, KC Steiner (2019). Adaptation to climate in five eastern North America broadleaf deciduous species: Growth clines and evidence of the growth-cold tolerance trade-off. Perspect Plant Ecol 37:64‒72

Li B, CN Dewey (2011). RSEM: Accurate transcript quantification from RNA Seq data with or without a reference genome. BMC Bioinform 12; Article 323

Liu Y, ZZ Xin, J Song, XY Zhu, QN Liu, DZ Zhang, BP Tang, CL Zhou, LS Dai (2018). Transcriptome analysis reveals potential antioxidant defense mechanisms in Antheraea pernyi in response to zinc stress. J Agric Food Chem 66:8132‒8141

Livak KJ, TD Schmittgen (2001). Analysis of relative gene expression data using real-time quantitative PCR and the 2⁻^ΔΔCTmethod. Methods 25:402‒408

Ma Y, X Dai, Y Xu, W Luo, X Zheng, D Zeng, Y Pan, X Lin, H Liu, D Zhang, J Xiao, X Guo, S Xu, Y Niu, J Jin, H Zhang, X Xu, L Li, W Wang, Q Qian, S Ge, K Chong (2015). COLD1 confers chilling tolerance in Rice. Cell 160:1209‒1221

Mantyla E, V Lang, ET Palva (1995). Role of abscisic acid in drought-induced freezing tolerance cold acclimation and accumulation of LTI78 and RAB18 proteins in Arabidopsis thaliana. Plant Physiol 107:141‒148

Mao X, T Cai, JG Olyarchuk, L Wei (2005). Automated genome annotation and pathway identification using the KEGG orthology (KO) as a controlled vocabulary. Bioinformatics 21:3787‒3793

Paul MJ, SP Driscoll (1997). Sugar repression of photosynthesis: The role of carbohydrates in signaling nitrogen deficiency through source: Sink imbalance. Plant Cell Environ 20:110‒116

Perochon A, D Aldon, JP Galaud, B Ranty (2011). Calmodulin and calmodulin-like proteins in plant calcium signaling. Biochimie 93:2048‒2053

Rajkowitsch L, D Chen, S Stampfl, K Semrad, C Waldsich, O Mayer, MF Jantsch, R Konrat, U Blasi, RRN Schroeder (2007). A chaperones RNA annealers and RNA helicases. RNA Biol 4:118‒130

Rooy SSB, GH Salekdeh, M Ghabooli, M Gholami, R Karimi (2017). Cold-induced physiological and biochemical responses of three grapevine cultivars differing in cold tolerance. Acta Physiol Plantarum 39:264-276

Ruelland E, MN Vaultier, A Zachowski, V Hurry (2009). Cold signalling and cold acclimation in plants. In: Advances in Botany Research, Vol. 49, pp: 35‒150. Kader JC, M Delseny (Eds.). Elsevier Ltd.

Shinozaki K, K Yamaguchi-Shinozaki (2000). Molecular responses to dehydration and low temperature: Differences and cross-talk between two stress signaling pathways. Curr Opin Plant Biol 3:217‒223

Tan W, QW Meng, M Brestic, K Olsovska, XH Yang (2011). Photosynthesis is improved by exogenous calcium in heat-stressed tobacco plants. J Plant Physiol 168:2063‒2071

Tatusov RL, MY Galperin, DA Natale (2000). The COG database: A tool for genome scale analysis of protein functions and evolution. Nucl Acids Res 28:33‒36

Thomashow M (2001). So, what’s new in the field of plant cold acclimation? Lots! Plant Physiol 125:89‒93

Trapnell C, BA Williams, G Pertea, A Mortazavi, G Kwan, MJV Baren, SL Salzberg, H Upadhyaya, SK Panda, BK Dutta (2011). CaCl₂ improves post-drought recovery potential in Camellia sinensis (L.) OKuntze. Plant Cell Rep 30:495‒503

Wang WH, EM He, Y Guo, QX Tong, HL Zheng (2016). Chloroplast calcium and ROS signaling networks potentially facilitate the primed state for stomatal closure under multiple stresses. Environ Exp Bot 122:85‒93

Williams JP, MU Khan, K Mitchell, G Johnson (1988). The effect of temperature on the level and biosynthesis of unsaturated fatty acid in diacylglycerols of Brassica napus leaves. Plant Physiol 87:904‒910

Xie JL, JQ Qi, XY Huang, N Zhou, Y Hu (2015). Comparative analysis of modern and ancient buried Phoebe zhennan wood: Surface color chemical components infrared spectroscopy and essential oil composition. J For Res 26:501‒507

Yang TY, XS Huang (2018). Deep sequencing-based characterization of transcriptome of Pyrus ussuriensis in responses to cold stress. Gene 661:109‒118

Young MD, MJ Wakefield, GK Smyth, A Oshlack (2010). Gene ontology analysis for RNA-Seq: Accounting for selection bias. Genome Biol 11:14-25